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Abstract 

We construct efficient Monte Carlo updating algorithms for two classes of pure SU (N) lattice 
gauge actions with non-linear dependence on the link variables. Our construction generalises 
the method of auxiliary variables used by Fabricius and Haan in the framework of Eguchi- 
Kawai models. We first review the original Fabricius-Haan method of constructing a pseudo- 
heatbath algorithm for fully reduced models, and discuss its extension to lattices with any 
number of reduced directions. We then use a similar method to construct updating algorithms 
for generic SU(N) mixed Wilson actions. We construct explicit examples of algorithms for 
Wilson actions whose plaquettes are in an irreducible representation of SU (N) with iV-ality 
k < 3. We also construct updating algorithms for the lattice version of centre-stabilised 
SU(N) Yang-Mills theories defined on R d_1 x S 1 , including the case of a fully reduced compact 
direction. We simulate the new algorithms and show that they are, in general, significantly 
more efficient than their Metropolis counterparts. 
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1 Introduction 

Simulations of pure lattice gauge theories normally involve the generation of Markov chains of 
gauge field configurations with the help of a local Monte Carlo algorithm. Such configurations 
(link variables on hypercubic lattices) are generated with respect to a probability distribution 
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characterised by the Boltzmann weight of the Euclidean partition function. Configurations are 
generated locally, i.e. the link variables are updated individually while keeping the remaining 
links fixed. A new configuration results from a sequence of local updates over the whole lattice 
(one sweep). It is thus necessary to know the form of the Boltzmann probability distribution 
restricted to only one link. The purpose of the local Monte Carlo (MC) algorithm is to sample 
this distribution function efficiently However, save in some special cases, the direct sampling of 
the probability distributions of individual links is often impracticable. Therefore, it is essential 
to find a good algorithm that samples such distributions efficiently, albeit indirectly. The ideal 
MC algorithm should be ergodic, fast, and produce decorrelated data. 

One straightforward possibility is to generate new link variables using a local Metropolis 
algorithm p] . Metropolis has the advantages of having an universal scope (it can sample any 
probability distribution) and of being very easy to implement. However, it is usually slow, 
decorrelates poorly, and may fail to be ergodic. Also, it requires the acceptance rates of new 
link proposals to be tuned at optimal values (which must be at around 50%, in order to avoid 
'too fast' or 'too slow' an exploration of the configuration space). Therefore, despite the ad- 
vantages, Metropolis is not a very efficient algorithm for simulating lattice gauge theories. In 
many situations, however, it is the only known algorithm. Fortunately, for pure SU(N) lattice 
gauge theories with the standard Wilson action it is possible to construct a more efficient 
alternative to the Metropolis algorithm. 

Consider the standard SU(N) Wilson action (up to an additive constant) with plaquettes 
in the fundamental representation: 

S F (f3 F ;[U\) = -^5>eTr{[/ p } (1) 

p 

Here p labels positively oriented plaquettes in the <i-dimensional hypercubic lattice A, U p is 
the plaquette operator, defined as the path-ordered product of all link variables U^ x G SU(N) 
in the boundary of p, 

U p ee U^U^U^Ul (2) 

(/t denotes the unit lattice vector in the /^-direction), f3 F is the bare lattice coupling, and Tr de- 
notes the character of the fundamental representation of SU(N). The probability distribution 
of link configurations [U] is Boltzmann: 

dP(/3 F ;[U}) = fi n [U] exp [U])) /Z F (J3 F ) (3) 

where /Xh denotes the product of the SU (A) -invariant Haar measures of all link variables, 

d 

mm = nn^ (4) 

and Zp is the Euclidean partition function of the lattice gauge theory: 

Z f {Pf) = [MU] exp (S F (P F ; [U])) (5) 
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Because the Wilson action Q is linear with respect to each link variable, the probability 
distribution of individual links is given by: 

dP{U„ tX ) oc dU^ x ex V (ReTr{Vl x U^ x }) (6) 

where V^ x is the sum of all 'staples' that couple to multiplied by the lattice coupling: 



Pf 



d 



u = \ 

It is important to note that V^ jX does not (and should not) depend on U^ x . 

In the SU(2) case, it is possible to perform a very efficient sampling of the probability 
distribution ^ |21 El S]- It relies on the fact that any sum v of SU{2) matrices, e.g. ([7]), is 
always proportional to an SU (2) matrix. Because of this property, the probability distribution 
of individual SU(2) matrices is of the form: 



dP(u) oc du exp ReTr j (t>£ 



where u,v^ x G SU(2), £ = 1 /det(w), and du is the S , f/(2)-invariant Haar measure. The 
S77(2)-projected sum of 'staples' v^ 1 can be absorbed by the Haar measure. Consequently, 
the probability distribution d8l) only depends on ao = Tr {u} G [—1, 1]: 



dP(ao) oc da (l — Oq) 2 ex P i^ a o) (9) 

which can be sampled very efficiently. Since SU(2) is isomorphic to S* 3 , ao corresponds to a 
'latitude' parameter. The traceless part of u is then uniformly distributed over the S 2 'equator', 
i.e. with a probability distribution given by the solid angle d 2 Q, which can be sampled trivially. 
The algorithms [2JE1I1] that generate random SU(2) matrices with the probability distribution 
(|8| are the heatbath algorithms. 

There are two other very useful SU(2) algorithms for lattice gauge theory simulations: 
the overrelaxation and cooling algorithms. The overrelaxation algorithm j5] performs 'large' 
local changes to the gauge field configurations that keep the Boltzmann weight invariant. 
Because the sum v of SU(2) matrices is always proportional to an SU(2) matrix, an exact 
overrelaxation update of u is obtained with the proposal: 

u i — y u = (vC 1 ) ■ u ] ■ (vC 1 ) e SU{2) (10) 
Using overrelaxation updates together with heatbath updates improves the efficiency of the 



overall MC algorithm, because (10) corresponds to take a large 'step' on SU{2) as compared 
with the typical heatbath 'steps'. This allows a broader exploration of the configuration space 
around the local minimum probed by the simulation. On the other hand, the cooling algorithm 
|6] generates gauge field configurations that minimise the lattice action locally. Because the 
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sum v of 577(2) matrices is always proportional to an SU(2) matrix, an exact cooling update 
of u is obtained with the proposal: 

u i — y v! = (vC 1 ) e SU(2) (11) 

The cooling of lattice gauge fields suppresses their quantum fluctuations. This is useful for the 
search of topological structures in the vacuum of the gauge theory, e.g. instantons. 

For gauge groups larger than 577(2), however, the situation is different: the sum of SU(N) 
matrices is not proportional to an SU (N) matrix, in general. Therefore, none of the algorithms 
described above can be directly extrapolated to iV > 2. A way to circumvent this obstacle is 
to construct algorithms that only update the SU(2) subgroups of SU(N), and not the whole 
group [?]. Due to the lattice action ^ being linear with respect to individual links U^ jX , it 
is then possible to single out each of its iV(iV — l)/2 SU(2) subgroup elements u C U^x, 
and obtain a distribution for u that is of the form (J8|. In this way, the SU(2) algorithms 
described above can still be used. The update of all (or a subset^ of all) SU(2) submatrices of 
a link variable results in a new link that is compatible with the probability distribution ([6]). 
The algorithm that generates random SU(N) matrices is known as Cabibbo-Marinari pseudo- 
heatbath algorithm [7]. Even though it doesn't sample ^ directly, it performs much better 
than the Metropolis algorithm in all aspects: it provides faster equilibration times, smaller 
autocorrelations, and it doesn't require any tuning. Because of all this, it has naturally become 
the standard algorithm for the generation of equilibrium configurations in pure lattice gauge 
theoriesJ3 

In the same way, the overrelaxation and cooling of SU(N) gauge configurations can be 
done by restricting such updates to the 577(2) subgroups. Recently, algorithms for the overre- 
laxation of the full SU(N) group, and not just of its SU(2) subgroups, have been suggested |9J. 
They are based on the singular value decomposition of V^x, and are shown to perform better 
than the overrelaxation of £77(2) subgroups. In the same line of thought, a cooling algorithm 
for the full SU(N) group can easily be constructed using the singular value decomposition of 

However, the SU(N) updating algorithms discussed above only apply if the probability 
distribution of individual link variables is of the form (|6]), i.e. if the lattice action is linear 
with respect to each link. If not, the contribution to the partition function of individual 577(2) 

1 The subset of SU(2) subgroups must be such that the remaining subgroups do not generate a left ideal 
[TJ. The minimal number of SU(2) subgroups that must be updated is therefore N — 1, i.e. those of the form 

(uij represents the 2x2 submatrix whose diagonal elements lie on the positions i and j of the diagonal 
of the N x N matrix it belongs to) . Despite the large difference between the minimal number and the total 
number of SU (2) subgroups, which becomes very significant for large N, it is recommended to update all of 
them. Updating too small a number of SU(2) subgroups may result in a poor performance. In our simulations, 
we always update all the N{N - l)/2 ST/ (2) subgroups of SU{N). 

2 An exact heatbath algorithm has been suggested for ST/(3) by Pietarinen [SJ. However, it is hard to 
implement, and so the Cabibbo-Marinari pseudo-heatbath strategy has also been adopted as the standard 
heatbath algorithm for SU(3). 
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subgroups is not of the form ([8]), and so the SU{2) algorithms discussed above cannot be used. 
For the same reasons, the full SU(N) overrelaxation/cooling algorithms cannot be used, too. 
In sum, if a lattice action is not linear with respect to the link variables, and no other efficient 
algorithms are known, then using a Metropolis algorithm is the only way to simulate such a 
theory. 

In this paper, we construct MC updating algorithms for some SU (N) lattice gauge actions 
whose dependence on the link variables is nonlinear. Our strategy is based on a generalisation 
of the Fabricius-Haan method of auxiliary variables [3j . We perform numerical simulations with 
the new algorithms and show that they are, in most cases, significantly more efficient then their 
Metropolis alternatives. In Section [2j we review the Fabricius-Haan method of constructing 
a pseudo-heatbath algorithm for the Eguchi-Kawai model. We then generalise their method, 
and apply it to the case of hypercubic lattices with any number of fully reduced directions. 
In Section [3j we apply similar methods to pure lattice gauge theories with a mixed Wilson 
action. We explicitly construct MC updating algorithms for Wilson actions whose plaquettes 
are in irreducible representations of SU(N) with iV-ality k < 3. We then study the numerical 
performance of these algorithms against Metropolis. In Section |4j we deal with the lattice 
regularisation of centre-stabilised SU (N) Yang-Mills theories defined on R 3 x S 1 . We explicitly 
construct MC an updating algorithm for these theories and study its numerical performance 
against Metropolis. We conclude our paper with a discussion on the advantages, efficiency and 
applicability of the new algorithms. We summarise in Appendix [X] all the algorithms proposed 
in this paper. 



2 Reduced lattices 

A situation where the probability distribution of individual links is not of the form ^ occurs 
already in pure lattice gauge theory with the standard Wilson action. If one or more lattice 
directions of a hypercubic lattice are reduced (i.e. have length = 1 in lattice units, for 
some direction fi), the Wilson action ([T]) becomes quadratic with respect to the link variables 
orthogonal to the reduced directions. Because of the this, gauge field configurations on reduced 
lattices would have to be generated with a Metropolis algorithm. However, Fabricius and Haan 
[3] were able to construct an efficient pseudo-heatbath algorithm for the case of a fully reduced 
(l d ) lattice, also known as Eguchi-Kawai (EK) model [lOjj^They used a method of auxiliary 
variables, which we review below. 

2.1 Fabricius-Haan method 

The EK model [10J is the original proposal of a matrix model for the large N limit of pure 
SU(N) lattice Yang-Mills theories. Initially, it was believed that their planar sectors would 
coincide, but the spontaneous breaking of an important symmetry eliminated that hope \V2\ . 

3 In fact, the pseudo-heatbath algorithm of Fabricius and Haan was originally constructed for l d lattices 
with twisted boundary conditions, also known as twisted Eguchi-Kawai (TEK) models |TT]. We omit the twists 
factors in the lattice action because they are not relevant to our discussion. 
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Since then, a number of alternative reduced models have been suggested [T2| [TT| IT3] . The 
goal is to find a zero- volume model of the full gauge theory in the large N limit, in which 
the set of global symmetries relevant for the large iV equivalence stay intact for all values of 
the coupling. This is hard to achieve, however, because those symmetries are sensitive to the 
volume. 

The action of the EK model is simply the fundamental Wilson action Q on a l d lattice: 

Sek(Pf;[u]) = -jfY, ReTr { u » u " u t u t} ( 12 ) 



Notice that the total absence of spacetime degrees of freedom in (|12j) makes the plaquette 
operators quadratic on the link variables. Therefore, the probability distribution of individual 
links cannot be put in the form of (|6]), because the sum of 'staples' (J7| would not be indepen- 
dent of Ufj,. Consequently, none of the efficient updating algorithms discussed in Section [T] can 
be applied directly, and the Metropolis algorithm seems to be the only alternative to simulate 
the EK model. 

Fabricius and Haan [3] circumvented this no-go by adding auxiliary degrees of freedom to 
the EK model. The new lattice variables are complex N x N matrices associated with the 
plaquettes, Q^ v = V// < v. They are given a free dynamics completely decoupled from 
the gauge field. In other words, the spurious degrees of freedom Q^ u are random matrices 
with the normal distribution, whose only effect is to multiply the partition function of the EK 
model by a constant factor: 

Zek(Pf) oc j n n [U] exp(-S EK ((3 F ; [U])) x J [dQ^dQ] exp (~ J^Tr [q\ v Q^]\ (13) 

* J^- ' 

constant 

/xh[^] is the product of S , f/(A^)-invariant Haar measures of all d link variables, and [dQ^dQ] 
is the product of standard flat measures for the auxiliary fields, 

d N 

[dQ ] dQ] = W\{ dRe(Q^) ab d\m{Q^ v ) ab (14) 

ji<v a, 6=1 

To simplify the notation, we use \i G to denote the Gaussian measure with variance a 2 of a 
complex variable z G C, 

fi a (z) = dz*dz (27T(T 2 ) _1 exp ( \z\ 2 J (15) 

or more generally, of a complex N x N matrix A G M(N, C), 



fX a {A) = dtfdA (27ra 2 )- Jv2 exp (--^Tr{AU} 



(16) 
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In this notation, the Gaussian integral multiplying the EK partition function is: 



f MQ] = i (ir) 

where m[Q] = Ylt<v MQ»»)- 

The second step in the Fabricius-Haan method consists in performing a particular change 
of variables (Q, U) i— >■ (Q, U) given by: 

Qfj.v = f^Y {Q^-U^Uv-UvUj (18) 



This transformation keeps the integration measure (14) invariant, up to a constant factor. 
It also cancels out the quadratic terms in the EK action, replacing them with Gaussian and 
linear terms. To see this, consider the effect of this change of variables on the exponent of the 



Gaussian term in (13): 



Gaussian -)■ - ^ ^ Tr { Q^Q^ } 

\i<V 



linear + ^7 £ ReTr WP^} 



d 



N 



cancels out + ^Y1 ReTl i U » U v U l U l } ( 19 ) 

pL<U 



The change of variables produces Gaussian terms for the Quv, linear terms on the link variables 



and quadratic terms that have the same functional form as the terms in ( 12 ), but with opposite 



sign. Therefore, the quadratic terms in (13) are cancelled out and only the Gaussian and linear 



terms survive. The partition function of the EK model with auxiliary variables then becomes: 
Zek(Pf) = JfJm[U]fi*[Q] ex P (^J2 Re ^{QU U » U "}) ( 2 °) 



with cr 2 = N/[3p. Since the exponent of the Boltzmann factor in (20) is linear with respect to 
each link variable, the Cabibbo-Marinari pseudo-heatbath can be used to update them. The 
update is performed with respect to the probability distribution of individual links (|6]), where 
U^x = Up, and V^ jX = is analogous to the sum of 'staples' ([7]): 

d 

v, = J2 {utQi»> + Qi»>ut) ( 21 ) 
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Note that does not depend on U^. In essence, the change of variables in (18) is an example 
of a Hubbard-Stratonovich transformation |14] . 

Despite the apparent coupling between the gauge and auxiliary degrees of freedom in 



(20), the gauge dynamics is completely unaffected by the presence of the auxiliary fields: the 



partition function (20) can always be transformed back to its original form (13) via the inverse 



of the non-singular transformations (18). Another apparent paradox resides in the fact that we 
are actually increasing the number if degrees of freedom that need to be updated in numerical 
simulations. However, all components (Q^ab °f the auxiliary fields are independent normally- 
distributed complex numbers. These can be generated very fast with a known algorithm, like 
the Box-Muller transform^] Therefore, the time needed to update the auxiliary variables is 
negligible as compared with the time needed to update the link variables. Consequently, the 
Fabricius-Haan method for the EK model is actually faster than its Metropolis alternative, 
given in [16J. 

2.2 Partially reduced lattices 

The Fabricius-Haan method for the EK model can easily be generalized to lattices with any 
number of reduced directions]^] In such lattices, only those plaquettes that are parallel to 
at least one reduced direction are quadratic with respect to the link variables. We call them 
reduced plaquettes. The remaining unreduced plaquettes are linear, hence the respective action 
terms do not need to be replaced. It then suffices to introduce one normally-distributed N x N 
complex matrix per reduced plaquette, i.e. Q^v^x = Qu/i,x such that = 1 or N v = 1. The 
partition function ^ of the partially reduced lattice gauge theory is then multiplied by the 
Gaussian integral over those auxiliary variables, J fi\[Q], which again amounts to multiply the 
partition function by '1'. With the following change of variables: 

i 

Q[IV,X ^ jy- J {Q^LU.X Uft^xUi/^+fl Ui/^xU^^x+u) (23) 

all action terms involving reduced plaquettes are eliminated, and the partition function ^ 
becomes: 

Z F ((3 F ) = J fj, n [U] fi ff [Q] exp{-S' F {/3 F ;[Q,U])) (24) 



4 The Box-Muller transform 1 15] is a method for generating pairs of independent random real numbers 
(171,32) with normal distribution dP(gl 7 <?2) = dg\dgi (27r<7 2 ) -1 exp(— (gf + gf)/2) from a pair of random 
real numbers (^1,^2) uniformly distributed over the unit interval (0,1]. The pair (171,92) can also be un- 
derstood as a random complex number z = g\ + ig^ with normal distribution over the complex plane 
dP(z) = dz*dz (2ira 2 )~ 1 exp(— |z| 2 /2). The Gaussian complex number z is generated from the pair (u\,u%) 
via the Box-Muller transform: 

z = (-21nzti)5 exp(i27ru 2 ) (22) 

5 An example of the Fabricius-Haan method applied to partially reduced lattices is used in [17] . 
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where a 2 = N//3p, and S' F is the 'linearised' action: 



S' f (Pf;[Q,U}) = -§E E ReTr {tW 

(L, 1 ,L I/ >1) 

-^E E R^WWV^A + M^)} (25) 

(L M = 1VL„ = 1) 

The contribution from the unreduced plaquettes is unchanged, because those terms are 
already linear in the link variables. On the other hand, the quadratic contributions coming 
from the reduced plaquettes are replaced by the linear terms involving auxiliary variables. 
The Cabibbo-Marinari pseudo-heatbath can then be used to update the link variables. The 
updates are performed with respect to the probability distribution ^ of individual links, 
where the sum of 'staples' V^ x is given by: 

a d 



+ % E {Q^ U U» + Ut, x Q»v )X ) (26) 



i/=i 

d 



i/=l 



The Fabricius-Haan algorithm for partially reduced lattices is summarised in the Appendix 
|A.1| The Fabricius-Haan algorithm for the EK model is contained as the special case for which 
the lattice is fully reduced. 



3 Mixed actions 

A suitable action for a SU(N) lattice gauge theory must be gauge-invariant and converge to 
the SU (N) Yang-Mills action in the continuum limit. It is well known |18j that class functions^] 
on SU(N) can be used to construct lattice actions that satisfy the conditions above. Examples 
of class functions are the characters of the irreducible representations TZ of SU(N). These are 
used to define an important class of lattice actions, namely Wilson actions whose plaquettes 
are in different representations of SU(N): 

SniPn;[U]) = ~^E Re ^™ ( 27 ) 

6 A class function is a function / defined on a group G that is constant on the conjugacy classes of 
G, i.e. f{ghg~ 1 ) = f(h),\/g,h G G. Here, g corresponds to a gauge transformation at a lattice site, and h 
corresponds to a lattice operator. Therefore, the definition of class function corresponds to the statement that 
/ is gauge-invariant. 
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where xn is the character of the representation 1Z, dn = d-ji = XtzW is its dimension, and 
/3ft is its associated bare lattice coupling. Equally suitable are the lattice actions consisting of 
arbitrary linear combinations of irreducible characters, known as mixed Wilson actions: 

S mi M[U]) = ^MMU]) (28) 
n 

— # 

where (3 denotes the set of independent lattice couplings 

All irreducible characters of SU (N) can be expressed in terms of the character of the funda- 
mental representation, \f = Tr. In this paper we only consider explicitly those representations 



of SU(N) with iV-ality k < 3, whose characters are given by: 

XF (U) = Tr{U} (29) 

XA (U) = |Tr{f/}| 2 -l (30) 

X(2)(U) = 1{Tt{U} 2 + Tt{U"}) (31) 

X(i,i)(U) = ^(Tr{f/} 2 -Tr{f/ 2 }) (32) 

X(3)(U) = ^(Tr{[/} 3 + 3Tr{f/}Tr{t/ 2 }+2Tr{f/ 3 }) (33) 

X(i,i,i)(U) = ^(Tr{f/} 3 -3Tr{f/}Tr{f/ 2 } + 2Tr{f/ 3 }) (34) 

X(2,i)(U) = ^(Tr{t/} 3 -Tr{f/ 3 }) (35) 



Here U is a generic element of SU(N), and the label in the characters denotes the Young 
tableau of a representation (e.g. (2) = m, (1,1) = 0, etc.); F denotes the fundamental 
representation, and A denotes the adjoint representation. 

It is clear that only the fundamental Wilson action is linear with respect to link variables. 
Consequently, the efficient updating algorithms discussed in Section [T] cannot be used directly 
in different representations. For this reason, all numerical studies of mixed Wilson actions have 
been performed using algorithms based on Metropolis, like the Cabibbo-Marinari-Metropolis 
and the overrelaxation- Metropolis algorithms |19| . or the biased Metropolis algorithm for the 
SU (2) fundamental /adjoint action j2D] - However, it is possible to imagine that a generalisation 
of the Fabricius-Haan method could also be applied to mixed Wilson actions, and in that way 
evade the nonlinear terms coming from the SU(N) characters. 

The philosophy behind the Fabricius-Haan method consists in adding a minimal (but 
sufficient) number of normally-distributed auxiliary fields, together with an appropriate change 
of variables, in order to eliminate all the nonlinear terms in the lattice action. These are then 
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replaced by linear terms on the link variables and Gaussian terms for the auxiliary fields. The 
possibility of using the efficient updating algorithms of Section [T]follows immediately. We argue 
that, by choosing a fair number of auxiliary variables and the correct transformation rules, it 
is ideally possible to eliminate all nonlinear plaquette terms in arbitrary mixed Wilson actions 
(and eventually in any lattice action whose dependence on the plaquettes is polynomial). 
Some of the auxiliary fields may eventually generate new nonlinear terms that also need to 
be eliminated. But as long as these are 'less nonlinear' than the terms they replace (e.g. by 
depending on smaller powers of the plaquette), the total amount of 'nonlinearity' is reduced, 
the process eventually stops, and all nonlinear terms are eliminated after a sufficient number of 
auxiliary fields is introduced. The objective is then to keep the number of necessary auxiliary 
fields at a minimum. We will make these statements more precise below. 



3.1 Linearisation of arbitrary characters 

For each representation 1Z of SU(N), we wish to linearise the corresponding Wilson action 



(27) with respect to the plaquette operator, U p . By 'linearising' the lattice action we mean 
replacing its nonlinear terms with linear and Gaussian terms, after adding a sufficient number 
Tin of normally-distributed auxiliary variables. Effectively, this linearisation corresponds to a 
Hubbard-Stratonovich transformation of the corresponding partition function. 

Let us associate to each positively oriented plaquette on the lattice n-ji complex N x N 
matrix variables Qpl )X (i — 1, • • ■ , Utz, and fi < v) with the normal distribution Hi(Q^l,x)- For 
the negatively oriented plaquettes we define Qvl,x = Q^x- First, we multiply the partition 
function of the mixed lattice gauge theory by j /J>i[Q] = 1, where fii[Q] is the product of 
Gaussian measures of all auxiliary variables. We then make a change of variables of the form: 

Qf = y/2fc{Q® - h$>) (36) 

Here p = (nv,x) labels plaquettes, [3' n = finN a jd-ji is a redefinition of the lattice coupling]^] 
Qp are the transformed auxiliary variables, and h p ^ = h^(Q p ,U P ) are functions of the 
plaquette U p and of a single auxiliary variable Q p with j < i. The condition j < i, whose 
origin becomes clear below, ensures that the change of variables (Q,U) h-> (Q, U) is non- 
singular, hence invertible. In fact, its Jacobian is just a non-zero constant. In other words, 
the integration measure of the HS-transformed partition function is invariant under such a 



change of variables, up to a multiplicative constant. After the change of variables (36), all 



nonlinear terms coming from the Gaussian measure fix[Q] must have been exactly cancelled 



by the nonlinear terms in the Wilson action (27). In the end, only the linear terms on U p and 
the Gaussian terms for Q p must survive. 

The way we choose suitable h's is the most straightforward possible. Basically, we associate 
one auxiliary variable Q p ^ to each nonlinear term of xn- in order to eliminate them, we perform 
a change of variables Q p — > Q p using a h p that only depends on U p , i.e. hp^ = W\U P ). If the 



7 a is an integer often equalling the degree of d-ji as a polynomial in N minus one. See Table [l] for the 
definitions of for each particular representation. 
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h's themselves generate new nonlinear terms, these must be eliminated by adding even more 
auxiliary variables (again, one auxiliary variable per nonlinear term). Each of these secondary 
nonlinear terms depends both on U p and on the auxiliary variable Q p ^ that generated it. 
Hence the new h's that we need to eliminate them must also depend on U p and on Qp\ 
i.e. hp = h^(Qp ,U P ) with j < i. If the secondary nonlinear terms have a 'lower iV-ality' 
than the terms they replace (i.e. if the original and secondary nonlinear terms, respectively 
Oi(U p ) and 2 (U p ), transforms under U p -> zU p as z kl Oi(U p ) and z k2 2 (U p ), with z £ Z N 
and ki > k 2 ), the process may be repeated until all nonlinear terms are eliminated and a 
linearised Wilson action emerges. 

Let us split foW into its linear and nonlinear parts: 

h®(Q®,U p ) = h®{Q<p,U p ) + h±(U p ) (37) 

The linear piece h p \ in its most general form, is given by: 

h®(Q®,U p ) = A 1 U p + A 2 Ul + A 3 Tr{A,U p } + A 5 Ti{A 6 Ul} (38) 

where the A t = Ai(Q p ) are complex N x N matrices depending solely on one Q p \ with j < i. 



The change of variables (36) has the following effect on the Gaussian measures yUi(Qp' ) ) 



1 Tr (QWtgW) @ _^ T r {qMqW} + 2/^ReTr {Q^hf} - &Tr {hfhf} 



2 I 



-/3^Tr {Qf^Qf} + 2/^ReTr {Q p i)r hf} + S nhn (U P ) (39) 



In the expression above, S n u n collects all the nonlinear terms generated by h®. We hypothesise 
that ^niin is either the symmetric of the sum of some nonlinear terms of \Tii or they are 
secondary nonlinear terms that can be cancelled out with the addition of more auxiliary 
variables. In other words, we assume that a linearisation of the Wilson action is possible. 



Given (38), the linear term on the r.h.s. of (39) can be rearranged as follows: 

2^ReTr{gWt^)} ® 2/^ReTr {gf%} (40) 

where g p are complex N x N matrices that only depend on the auxiliary fields: 

g®(Q<P) = A\Qf+Q^A 2 + AlTr{AlQf} + A 6 Tr{A 5 Q p ^} (41) 
and Ai = Ai(Q p ^), j < i. After all n-ji auxiliary variables are added and the change of variables 



(36) is performed, the exponent in the Gaussian measure fi\[Q] becomes: 
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i=l p 

Gaussian ->■ - P'nJ^Yl ^{Qp^Qf > 

i=l p 

linear -> + 2/3^ ^ ReTrf/^t/p} 

p 

cancels out ->■ - V, R e Xrc {^p} (42) 

P 

where is a complex N x N matrix that only depends on the auxiliary variables: 

ff = Y,9? (43) 
i=i 

In sum, encodes the nonlinear dependence of xn on U„. 



Once the Wilson action (27) is linearised, the efficient updating algorithms discussed in 
Section [T] can be applied directly. The updates are performed with respect to the probability 
distribution of individual links where V^ x = V™ is the sum of 'staples': 

d 

^fi,x = 2/?7£ ^ ] (fiiVjxUvfxUfax+pUl^p + U\ x _pf u ^ x _ j> U^ :X ^yU u ^ x -i,j r j^ (44) 



i/=i 



and = /^Jj., V/i < v. It must be noted that the order in which appears in the 'staples' 
depends on the particular choice for the initial point of the plaquette operator (|2|. Changing 
the initial point of the plaquette to any other vertex would simply result in a relocation of 

within the 'staples'. For the particular case when 1Z = F we have / F = 1, so we have 
recovered the original algorithm. 

In the next Sections we linearise explicitly the Wilson action for each representation of 
SU(N) with iV-ality k < 3; Table [j] summarises the changes of variables necessary for each 
linearisation. We then discuss the case of mixed actions, with and without reduced directions. 
In the end, we perform numerical tests on some of these algorithms against Metropolis, in 
order to check their accuracy and performance. 



3.2 Adjoint representation 

The best studied example of a mixed Wilson action is the one involving both fundamental 
and adjoint plaquette terms, 

S f+a (Pf,Pa;[U]) = S F (P F ;[U]) + S A ((3 A ;[U]) (45) 
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The parameter space of (Pf,Pa) lattice couplings is known as the fundamental/adjoint plane. 
This action is important, for example, in the study of the role of centre of the gauge group 
in colour confinement, because the adjoint character Xa{U p ) is invariant under centre shifts, 
U p i — ^ zU p , z e Z N . 



Due to the quadratic nature of \A (30), it is not possible to use the efficient updating 
algorithms of the fundamental Wilson action, discussed in Section [1} A Metropolis algorithm 
is often used instead. In particular, combining a Cabibbo-Marinari-Metropolis algorithm with 
overrelaxation-Metropolis sweeps results in a rather good performance, at least for SU(3) |19j . 
Recently, a biased Metropolis algorithm with heatbath efficiency has been constructed for the 
SU(2) case |20j. In this Section we construct a pseudo- heatbath algorithm for the SU(N) 
adjoint Wilson action, 

Sa(PaM) = -J^Y,p\^ U v}\ 2 (46) 

for all N, using the method of auxiliary variables. But because the double-trace term is always 
positive, the cases for positive and negative adjoint coupling Pa must be considered separately. 



3.2.1 A > 



When /3a is positive, the adjoint Wilson action (46) is always negative. In order to cancel 



out the double-trace terms, it suffices to introduce one complex number per plaquette z p with 
the normal distribution fii(z p ). We then multiply the partition function of the adjoint Wilson 
theory Z a (Pa) by the Gaussian integral f fJ>i\z\ — 1. Consider the change of variables: 



.1 \ z p-^ TT i u p} 



(47) 



where P' A = P A N 2 /(N 2 



- 1). The effect of (47) on the Gaussian exponent of fii[z\ is: 

_ I V* 17 l 2 © 

2^p 1p1 

~$a y, \ z p\ 



Gaussian — >■ 



linear — y 



cancels out — > 



Pa 
N 2 



The last term in the r.h.s. of (48) cancels out all the adjoint terms of (46) exactly, as long 



as Pa is positive. It comes from conjugating hp = Tr{U p }/N with itself. If Pa < 0, then it is 
not possible to do the same with ^-variables only. We discuss this case later. The partition 
function of the adjoint Wilson theory then becomes: 



Za(Pa) 



(49) 
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where fi a [z] = Y1 p ^( z p) an d °" 2 = ^/P'a- Since the exponent of the Boltzmann factor is 
now linear, the link variables can be updated using the efficient algorithms of Section [TJ The 
updates are performed with respect to the probability distribution of individual links ([6]), 
where V^ >x = V^ x is given by (44) with = z p /N. And because z p and are just complex 
numbers, the cost of updating the auxiliary variables in this case is negligible as compared 
with link updates. 

3.2.2 A < 



In the case of negative (3a, the adjoint Wilson action (46) is positive. Therefore, the double- 



trace terms cannot be eliminated using (47). However, it is possible to linearise the action if 



we introduce different auxiliary variables, namely one complex N x N matrix per plaquette Q p 
with the normal distribution fii(Q p ). As usual, we multiply the partition function 
with the Gaussian integral J fi\[Q] — 1. Then we consider the change of variables: 



Q p = \[W A Q 



U p --Tt{U p }l 



(50) 



where (3' A = \(3a\N / {N 2 — 1). The effect of (50) on the Gaussian exponent of fii[Q] is: 

Gaussian ->■ ~P'a^2 ^ {QpQp} 

+2/3^ReTr j (q p - ^Tr {Q p } l) ^ U p 



linear — > 



cancels out 



(51) 



The last term in the r.h.s. of (51 ) carries the correct sign to cancel out the double-trace terms 



with negative (3a from (46). The partition function of this theory then becomes: 

Z A {(3 A ) = j MU] H*[Q] exp ^^ReTr | (q p - ^Tr {Q p } U p ) ) (52) 

where fi a [Q] = Tlp^criQp), arid o~ 2 = 1/(3' A - Individual links are then updated with respect to 
the probability distribution (J6|, with V^ >x = V px given by (44), and f p is the matrix factor: 



Q p --Ti{Q p }t 



(53) 



Since Q p and f p are N x N matrices, and not just complex numbers, it is expected for the 
(3a < case to be less efficient than the (3a > case. 
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3.3 Higher N— ality representations 



The case of higher iV-ality representations is more complicated, because there are more non- 
linear terms to be eliminated, and these depend on higher powers of the link variables. This 
means that a larger number of auxiliary variables is needed in order to eliminate all nonlinear 
terms. In the following, we explicitly linearise the Wilson actions for two- and three-index 
irreducible representations of SU(N) using the method of auxiliary variables. 



3.3.1 TZ = 2± 

Here we consider Wilson actions with plaquettes in the symmetric (+) or antisymmetric (— ) 
two-index representations of SU(N), respectively 1Z = (2) or TZ = (1, 1) : 

S ± (P ± ;[U]) = - N{ ^ ± 1} E Re ( Tr ™' ± Tr KD (54) 

We consider both situations of positive or negative lattice coupling /3±, whose sign we denote 
by a = sgn(/3±). One way of linearising (54) requires the addition of three auxiliary complex 
matrix variables per plaquette Q p \ i = 1,2,3, and the corresponding change of variables: 

= VW ± {Q?-\{u p + ^{ul}t)^ (55) 
Qf = VW±(Qf ] -\(u p ±^uij} (56) 
Qf ] = VW ± {Qf ] -\{u p -^{u p }ty^ (57) 

where f3'± = 2\j3±\/(N±l). The variable Q p 1 ^ is responsible for the elimination of the quadratic 
term Tr{U p } 2 in the action, and its replacement by a linear term. However, it also generates 

a secondary nonlinear term of the form |Tr{L/p}| 2 . The variable Q p 2 ^ eliminates the quadratic 
term Ti{U p } in the action and replaces it with a linear term, with no other side effects. Finally, 

the variable Q p ^ eliminates the term |Tr{?7 p }| 2 generated by Q^K This is achieved using the 
transformation (57), which is similar to the one used in the elimination of negative-coupling 



adjoint plaquettes (50). The remaining linear terms contribute to the sum of 'staples' (44) 



with the matrix factor: 



ft = \ (q? + ^Tr{Q«t} + Q( 2 ) ± ^ 2)t + gi 3) -^Tr{Q( 3 )}l) (58) 

However, the above choice is not unique. Is is also possible to linearise the Wilson action 
(54) by adding only two auxiliary complex matrix variables per plaquette. The corresponding 
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change of variables is: 



§« = Vm (q p 1] - \ (*Ul ± ±U P + i T r{f/ p }l^ (59) 

The variable Q p 1 ^ is responsible for the simultaneous elimination of both nonlinear terms from 
X±- However, it also generates a |Tr{?7 p }| 2 term, which is promptly eliminated by variable 
Q { p>. After some algebra, we find the contribution of this change of variables to the sum of 



'staples' (44) to be: 



St = \ ± ^ + + (^r) ' («? - >W? }») 

Albeit different, both situations describe the same theories. However, it is clear that the 
second method, with only two auxiliary variables, is superior to the first method: the smaller 
the number of auxiliary variables, the more efficient we expect the corresponding MC algorithm 
to be. 

3.3.2 K = 3± 

Here we consider Wilson actions with plaquettes in the symmetric (+) or antisymmetric (— ) 
three-index representations of SU(N), respectively 1Z = (3) or 1Z = (1, 1, 1) : 

S±(P±;U) = ~ N{N2 ±* N - 2) E Re ( Tr ™ 3 ± 3Tr Tr i u p}+ 2Tr W}) ( 62 ) 
We again consider both situations of positive or negative lattice coupling 0±, whose sign we 



denote by a = sgn((3±). One way of linearising (62) requires the addition of seven auxiliary 
complex matrix variables per plaquette Qp\ i = 1, . . . , 7, and the corresponding change of 
variables: 



17 



Q« = y^^--^(U p Tr{U p } + aTr{Ul}t)^ (63) 
Qf = Vm^Q P 2) -^(U P TT{U p }±aU^ (64) 

Q? = V / 2^(^ 3) -^K + ^ P t )) (65) 
= - - ^Tr{f/ p }l) (66) 

Qi 5) = (gf - ig«C5 - ^Tr{f/ P }l) (67) 

Qf = \fW± (q p 6) - -^Qp )u I - Jj U p) (68) 

Qi 7) = Vm (q? ] - (§) 2 (tfp - ^Tr{[/ P }l^ (69) 

where ^ = 6|/5 ± |iV/(A^ 2 ± 3N + 2). The variable Q p l) eliminates the cubic term Tr{[/ p } 3 
from the action, but it also generates nonlinear terms of the form Th:{Qp Up\ r Tt{U p \ and 
|Tr{[/p}| 2 . The variable QP eliminates the cubic term Tt{U p }Tt{U p } and generates similar 
nonlinear terms, namely Tr{g p 2 ^?7p}Tr{[7p} and |Tr{[/p}| 2 . The variable Q p ^ eliminates the 
cubic term Tr{[/p} and generates one nonlinear term of the form Tr{Q p 3 ^U p }. The variables 
Qp 4 \ Qp an d Qp eliminate the secondary nonlinear terms generated by Q p X \ Q p 2 ^ and Q p 3 \ 
respectively, except the term of the form |Tr{£7p}| 2 ; in fact, the variables Q p ^ and Q p ^ also 

(7) 

produce such a term. The sum of these double-traces is eliminated by the variable Q p . All 
the seven auxiliary variables generate linear terms, whose contribution to the sum of 'staples' 



(44) is given by: 

f ± = ° Tr{0 (1)t |l ± — Q (2)t + ° C (3)t + — Tr|g (4) |l + — g (4)t c (1) 



x i2 y v p N' 

Just like in the previous case of two-index representations, the choice above is not unique. 
In fact, there is a cheaper way to eliminate all nonlinear terms, involving only five auxiliary 
variables. The corresponding change of variables is: 
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Q« = ^(Q^--^(U p Tr{U p } + aTr{Ul}t±3aU^ (71) 

= ym(Q?-^W + *U}j) (72) 
Q< 3 > = JW ± (of - ^Q { pH - jjTr{U p }l) (73) 

Qf = JW± [Qf - - ^p) (74) 

5f» = ^(^-( Z ^)*(^-^p}l)) (75) 

The variable Q p 1 ^ eliminates simultaneously both the Tr{U p } s and Tv{U p }Tt{U p } terms in the 
action, but it also generates nonlinear terms of the form Tr{Qp 1 ^[/p}Tr{C/ p } and |Tr{[/ p }| 2 . 
The variable Q p 2 ^ eliminates the cubic term Tr{?7p} and generates one nonlinear term of the 
form Tr{Q p 2 ^U p }. The variables Q^p and eliminate the nonlinear terms generated by 
QP and Q p 2 \ respectively, except the term of the form |Tr{t/ p } | 2 (which is also generated by 
Q P 3) ). The sum of these double-traces is eliminated by Q p 5 K In the end, the contribution of 



the remaining linear terms to the sum of 'staples' (44) is given by: 



+^ W + ^ 4 > + (^) 1 (flf - (™> 

3.3.3 ft =(2,1) 

Finally, we consider Wilson actions with plaquettes in the representation 1Z = (2, 1): 

s<m)(pm,v) = - N ^T_ 1} E Re ( Tr ™ 3 - Tr Wi) (77) 



The sign of /3(2,i) is denoted by a = sgn(/3±). One way of linearising (62 ) requires the addition of 



five auxiliary complex matrix variables per plaquette Qp \ i = 1, . . . , 5, and the corresponding 
change of variables: 
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Qp ] = v^^(^ 1) "7^(^ Tr ™ + aTr ^ t > 1 )) ( 78 ) 

Q? ] = y/^)(Q?-^W- aU i)) W 

Qf = (Q? - ^Q^U} - ^Tr{U p }l^ (80) 

Q\f ] = xj^^-^ftH-^) (si) 

Q ( p ] = J^(q® -^(u p -±Tr{U p }l)) ( 82 ) 



(2>1) V p V N 

where /3| 2 ^ = 3|/3(2,i)|iV/(iV 2 — 1). The variable eliminates the cubic term Ty{U p } 3 in the 
action, and it generates nonlinear terms of the form Tr{Q p U p }Tr{U p } and |Tr{f/ p }| 2 . The 
variable Q p 2 ^ eliminates the cubic term Tr{[/p} and generates one nonlinear term of the form 
Tr{Q p 2 ^Up}. The variables Q p 3 ^ and Q p ^ eliminate the secondary nonlinear terms generated 
by Qp 1 ^ and Q P 2 \ respectively, except the term of the form |Tr{[/p}| 2 (which is also generated 
by Q p °>). The sum of these quadratic terms is eliminated by the fifth variable. In the end, the 



contribution of the remaining linear terms to the sum of 'staples' (44) is given by: 



f ™ = i^ TrW ' 1,,}1 -i]? ? ,,+ ^ Tr{ ^ ,}1+ ^ <^ ™ <^ " , 

+1qm + ^ef'ef + ^ («« 5) - ^{Q?n) (83) 

Also in this case it is possible to find a smaller set of auxiliary variables that does the same 
job. In fact, only three auxiliary variables suffice to eliminate all the nonlinear terms from the 



action (77). The corresponding change of variables is: 



Qp ] = sJW^) (q p 1) ~ i (U P (Tr{U P }t ~ U p ) + a (U p + Tr{U p }lf)\ (84) 



1 1 



Q? = yj2P[ aA) \^P-^QWui--(U 9 -Tr{U p }l)j (85) 

AN -6\* / I 



= V 2/3 (^) K 3) " (r^rj /v' IV|/ 11 ) I ' s()) 

The variable Qp 1 ^ eliminates both nonlinear terms in the action, but it also generates non- 
linear terms of the form Tr{Q^U p }Tr{U p }, Tr{Q P 1)j U^} and \Tt{U p }\ 2 . The variable Q ( p ] 
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eliminates all secondary nonlinear terms generated by Q p : except the term |Tr{f/ p }| 2 , which 

(3) 

is also generated by it. The sum of these double-traces is eliminated by Q p . In the end, the 



contribution of the remaining linear terms to the sum of 'staples' (44) is given by: 



f(2,l) 
J V 



a 



Vqn 



(1) 



TrjQW}!) 1 



1 

N 



(Qf-TrlQfjl) 



— c (1) c (2) 



AN - 6\ 5 
3N 



N 



Tr{Q^}l 



3.4 Mixed Wilson actions 

The generalisation to the case of mixed Wilson actions is trivial, because a linearised mixed 
Wilson action is simply the sum of linearised Wilson actions for each representation: 



n 



Therefore, the probability distribution of individual links associated with a generic mixed 
Wilson action is given by ([HJ, where V^ x is the sum of 'staples': 

d 



i/=l 



and f^x = fp encodes the information about all representations involved: 

n 

The possibility of using the efficient MC algorithms discussed in Section [T] follows immediately. 



The MC algorithm for a general mixed Wilson action is summarised in Appendix A.2 

This straightforward approach to mixed Wilson actions may not be the most efficient, 
however. In this approach, each representation is treated independently, and the total number 
of auxiliary variables is n m i x = Yln 71 ^-- But in some situations it is possible to reduce n m i x . 
For example, consider the mixed L A + (2)' Wilson action. During the linearisation of X(2), 
a secondary nonlinear term of the form |Tr{[/ p }| 2 is generated, and it can be eliminated in 
simultaneous with the double-trace terms of Xa{U p ). Therefore, only two auxiliary variables 
per plaquette are needed, instead of the naive ua + n<2) = 3. This is the same ambiguity 
that exists in the case of individual representations, as shown above. In sum, the process of 
linearisation is not unique, and may be improved by a wise choice of auxiliary variables and 
respective change of variables. In the end, one must always choose the set with the smallest 
number of auxiliary variables. 
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3.5 A note on reduced mixed models 



A special note must be taken in the case of mixed Wilson actions on lattices with reduced 
directions. These may be useful to study the (non-)universality of the symmetry breaking 
transitions that invalidate the large TV equivalence in TEK models [2 lj . 

As usual, the construction of an efficient updating algorithm for a mixed reduced model 
consists in adding enough auxiliary variables in order to eliminate all the nonlinear terms 
in the action. In reduced models, however, the nonlinearities have two origins: the quadratic 
nature of the reduced plaquette operator, and the nonlinear nature of the SU(N) characters. 
Clearly, the nonlinearities associated with SU(N) characters need to be dealt with first. For 
that end, the linearisation of the mixed lattice action proceeds exactly as described in Sections 

mm 

The problem occurs in the final step, when dealing with the reduced plaquette operator. 
One would be tempted to solve the problem with the original Fabricius-Haan treatment for 
the EK model. The nonlinear terms to be eliminated are of the form: 

- 2/3^ ReTr {f^U^U^U^U^} (91) 

which differ from the EK plaquette terms by the matrix factor fj^ x . If one adds an auxiliary 
variable R^ v ^ x and perform the change of variables 

the linearisation with respect to the link variables is promptly achieved: 
--Tr/# R \ H 

Gaussian^- — /3' n Tl {R^ u ^R^^} 

non-Gaussian ~^Tr {f%}, x f™,x} 

linear +2/^ReTr {R^ x {J%.r ,, J . + U„ jX U^ x ) } 

cancels out -2/^ReTr {^U^U^U^Ul^} (93) 

However, new terms of the form Tr {fj^} x fj^, x } appear in the linearised action. Even though 
they do not depend on the link to be updated, they give a non-Gaussian weight to the auxiliary 
variables. Therefore, the efficient algorithms discussed in Section [T] cannot be applied, unless 
such terms are eliminated too. The elimination of the term Tr { f^xf^vx] * s no ^ eas Y- 111 
general, has a rather complicated dependence on the auxiliary variables, as can be seen 



in the examples derived in Sections |3.2f|3.4| The only situation for which we have an easy 
solution to the problem (or, at least, a solution that does not undermine the efficiency of the 
resulting algorithm) is the case of the adjoint representation. 
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For positive (3a, the solution is trivial. Given that f£, x = z^ UtX /N, the non-Gaussian term 
is in fact Gaussian: 

-^{f^XuJ = -flwl 2 (94) 

and so the problem is solved. The update of the link variables for a partially reduced adjoint 
Wilson action is performed with respect to the probability distribution of individual links ([6]), 
with V^ x = V£ x given by: 



V A _ 2 ^A_ ST^ ( 7 tj tj rrt I 7 * rrt tj tj 

i/=l 

+ 2 ^ E {jf^^M^n + U l R ^*) ( 95 ) 



The first term in the r.h.s. of the equation above is the contribution from the unreduced 
plaquettes, and the second term is the contribution from the reduced plaquettes. 

For negative (3a, f^ ux is a matrix. In order to eliminate this term, we first expand it in 
terms of the Q- variables: 



(5.3) 



Gaussian^ ~ P'a^ {Qfiv,xQ n^,x} 



ft 



nonlinear^ +^ |Tr {Q^W (96) 

This term is clearly non-Gaussian. However, it is the sum of a Gaussian term and a double-trace 
term that can easily be eliminated. For that end, we introduce yet another auxiliary variable 
^iiv,x with the normal distribution /i 1 (M Ml/ . r ), associated with positively-oriented reduced 
plaquettes (for negatively-oriented plaquettes we define M UflfX = M \ ux )- Let us perform the 
change of variables: 

M mx = V^^,x-^Tr{Q^ x }U^ x+0 Ul x+ ^ (97) 



The effect of (97) on the Gaussian exponent in /J,i(M fJ , l/iX ) generates the following terms: 
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Gaussian ->■ -P'a^ t { M U,x M ^x } 



linear +^ReTr {Tr {Q^} M^U^U^ x+ ~] 



fit 

cancels out ->■ --^ |Tr {Q^aJI (98) 

the last of which cancels out the nonlinear term coming from the non-Gaussian piece, thus 
solving the problem for negative /3 a- The sum of 'staples' for the probability distribution of 



individual links (44) is then given by: 



a 



+ 2/3^4 ^ ] {f^u,x^-^,xUl x+ ^ + Ul x Rp UtX + ^Tr {Q^^} M^xU^x+fa J 



(99) 



where f£ v x is defined as in ( 53 ) . 



For higher representations, the analogous of the expansion (96) of the non-Gaussian term 
results in a higher number of complicated nonlinear terms. In order to eliminate them, many 
more auxiliary terms would have to be introduced. This would certainly render the resulting 
algorithms inefficient and useless. Although unlikely, it is not a priori impossible that a clever 
choice of a smaller set of auxiliary variables would result in viable updating algorithms for 
higher-representation reduced models. 

3.6 Numerical tests 

We simulated some of the new algorithms proposed above, with the purpose of comparing them 
with the Metropolis algorithm. We considered mixed l F + TV Wilson actions in d — 4, i.e. 



the Wilson action (27) with plaquettes in the fundamental representation, F, and in another 



representation 1Z ^ F of SU(N). We simulated each theory with both a Metropolis algorithm 



and the new MC algorithm proposed in Section 3.4 We tested both algorithms for their 
compatibility, by checking if the expectation values of gauge-invariant observables coincide in 
both cases. We also tested them for their relative efficiency, by comparing the magnitude of 
the autocorrelations in the Markov chains they generate. 

For the thermal Metropolis updates we used a variant of the (1-hit) Cabibbo-Marinari- 
Metropolis (CMM) algorithm described in [19], and appropriately adapted to an arbitrary 
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representation TZ. In our CMM algorithm, new link proposals are generated via Cabibbo- 
Marinari updates with respect to the fundamental part of the action only, in which j3p is 
replaced with the tuning parameter of the Metropolis algorithm /3m (Avi < Pf)- The link 
proposal is then accepted or rejected a la Metropolis with respect to the full l F + TV action. 
The acceptance rates are tuned to stay in the range 40-60%. For the overrelaxation updates, 
we adapted the overrelaxation-Metropolis algorithm also described in [19]. In this algorithm, 
SU (2)- or SU (iV)-overrelaxed link variables are accepted or rejected with respect to the S-n 
part of the action only. Acceptance rates for the overrelaxation-Metropolis algorithm stayed 
well above 85%. 

For both the Metropolis and the new algorithms, each configuration update consisted of one 
thermal update followed by 5 overrelaxation updates. For each configuration, we evaluated the 
characters \f and \n of the plaquette, in order to estimate their expectation values. In each 
simulation we performed O(10 5 ) measurements, after discarding the initial 2,000 configurations 
for equilibration. The simulation parameters and measured observables, together with their 
naive confidence intervals, are shown in Table [2j 

The parameters of the simulations of the 'F + A action were chosen to coincide with 
those of |19j . In this way, we could compare our results with the literature. We considered two 
particular values for the adjoint lattice coupling, one positive and one negative, that were also 
considered in that article. We obtained compatible results for both plaquette characters, which 
is good evidence that our new algorithm reproduces accurately the l F+A* lattice gauge theory. 
For higher representations, we performed simulations of the l F + TV theories on a 8 4 lattice. 
In each of these cases, the observables calculated for a particular value of the lattice coupling 
also matched in both algorithms. There is, however, some discrepancies in the last significant 
digits of each observable. This is probably due to the fact that the confidence intervals shown 
in Table [2] do not take autocorrelations into account, i.e. they are underestimated. 

We also measured the autocorrelations of the Markov chains generated by both types of 
algorithms. We calculated the values of the fundamental trace of the plaquette m as a function 
of the MC time i, and then used them to estimate the normalised autocorrelation function: 

From C(t) we estimated the time-dependent integrated autocorrelation time: 



'hit 



(t) 



(101) 



i=i 



The estimator of the integrated autocorrelation time, T- m %, is the plateau value of Ti nt (t) when 
t — > oo. In terms of these autocorrelations, the new pseudo-heatbath algorithms performs 
significantly better than their Metropolis counterparts, as can be seen in FigsjT 12 The graphs 
show that the new algorithms decorrelate faster than Metropolis, for all representations, both 
in MC time and in effective CPU time. 
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We also observed that, except for the case of the adjoint representation, one pseudo- 
heatbath update takes longer to perform (in CPU time) than one Metropolis hit. This is 
not surprising, taking into account the large number of auxiliary variables that some repre- 
sentations require. Only for the adjoint case, which requires only one auxiliary variable per 
plaquette, were the pseudo-heatbath updates faster than Metropolis. In particular, this dif- 
ference was largest for 0a > 0, as expected. Nevertheless, the fast decorrelation of the new 
algorithms overshadows the disadvantage of slower updates, making them a superior alterna- 
tive to Metropolis. 

4 Double-trace deformations 

Recently, Unsal and Yaffe suggested |13| a strategy to prevent the spontaneous breaking of 
the centre symmetry in pure SU(N) Yang-Mills theories on manifolds with small compactified 
directions, namely M d ~ k x (S' 1 ) fe . The strategy consists in deforming the Yang-Mills action 
with a centre-stabilising potential, which contains double-trace operators dependent on the 
holonomies wrapping the compact directions. This mimics the effective potential contribution 
of adjoint fermions with periodic boundary conditions in the compact direction, which are 
known to stabilise the centre symmetry [22]. The gauge theories deformed by such double- 
trace terms are called centre-stabilised Yang-Mills (CYM) theories. 

The breaking of the centre symmetry is the reason why the EK model (and most of its 
variants) is not equivalent to pure Yang-Mills theories in the large N limit. CYM theories, 
however, are believed to be completely volume-independent in the large N limit. If this is 
true, then it would be possible to represent large N Yang-Mills theory defined on R d by 
the zero volume limit of the compact directions of CYM theories. Taking the volume of the 
compactified directions to zero with no consequences for the large N equivalence has clear 
analytical and numerical advantages. 

In this section we consider the lattice regularisation of CYM theories defined on manifolds 
with a single compactified direction, R, d_1 x S 1 . Due to the double-trace terms, the CYM 
lattice action is highly nonlinear with respect to the link variables, hence Metropolis seem 
to be the only possible algorithm to simulate such theories. However, we will show that it is 
possible to linearise the CYM lattice action using the method of auxiliary variables, and from 
there obtain a rather efficient MC algorithms. 

4.1 CYM theory on E d 1 x S 1 

The lattice action for the SU(N) CYM theory defined on a manifold with only one compact 
direction, namely R d_1 x S 1 , is given by: 

ScYM(a,p F ;[U}) = S F ((3 F ;[U]) + S dci (a;[U]) (102) 

where S F is the fundamental Wilson action ([I]) and a = («i, . . . , a|_Ar/2j) are free parameters 
of the CYM model that must be positive. The deformation potential Sdef is a sum of double- 
traces: 
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LJV/2J 



S dei (a;[U}) = JL_ E |Tr{^}| 2 (103) 

* xeAi n=l 

where A^ is the size of the compact direction (S 1 ) in lattice units, x labels the lattice sites 
on the (d — 1) -dimensional lattice Aj_ (the discretization of R d_1 ), and fi x is the holonomy 
(Polyakov loop) wrapping the compact direction, d: 

N t -1 

^ = -PU U ^= U d) *U dil[+r ..U diX+{Nt _ 1)S (104) 

i=0 



The CYM action (103) is highly nonlinear with respect to the link variables, especially for 
large A" (because it contains terms that depend up to the [A r /2j-th power of the link variable). 
In addition, Sd e f contains (^(A^ 1 ) terms, which makes the theory even harder to simulate in 
that limit. From these facts, we should expect to use a relatively large number of auxiliary 
variables in order to linearise the whole action. Let us consider the following set: 

Rn,x , 1 < n < K 
Q$ , 1 < m < n < K 
45 , 2 <n < K 

where K = |_7V/2j . We then multiply the partition function of the CYM theory by the Gaussian 
integrals J fi\[R, Q] = 1. We now show that the following change of variables linearises the 
CYM action: 



Rn,* = ( |S) 2 (Vx - (n» - ^Tr{Kn ) ) (105) 



QS = (^) 5 »& ) -»Sr 1) ni + ^)) ( 106 ) 
QSS = (fpr) 2 (« } - ^{RnM) (107) 

where we define Qi° x to be the traceless part of -R„ iX : 

Q { & = Rn,* - ^Tr{i? niX }l (108) 

The i?-variables cancel out the double-trace terms from the CYM action. We can see this 
when we apply the change of variables (105) on the exponent of the Gaussian measure fii[R\: 
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n=l 



N K 



Gaussian — > — d-1 

* n=l 

2iV 



linear — y 



T axReTr |gg f a 



2N K 



nonlinear +-^ T £ a n ReTr {Qgfnj} 

* n=2 
1 ^ 

cancels out -»■ +^TT E «n |Tr {^}| 2 (109) 



1 n=l 

The nonlinear terms in the r.h.s. of the expression above are cancelled by the Q-variables: 

K n-l 



n=2 m=l 



K n-l 

Gaussian E E «« Tr WS^K I 

* ?i=2 m=l 

7V * 

Gaussian -). --^ £ « n Tr {Q£ x ~ WgC^ D } 

t n=2 



non-Gaussian — y 



linear — > 



N K 



4—1 ^ ] tt n Tr {Qn^QnL.} 



2N K n ~ l 



i E E a - ReTr {QipQfc 1 ^} 



d 

n=2 m=l 



2N K 



linear -). £ « n ReTr {Q^Oj 



1 n=2 

2N - 

cancels out "T^T E «» ReTr (HO) 

4 n=2 

The only remaining pathological terms that need to be eliminated are non-Gaussian. This 
situation is similar to the case of reduced mixed models with a negative adjoint coupling (3a 
(see Section 3.5). Hence let us expand the non-Gaussian piece in terms of the /^-variables: 
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1 ^ ] {Qn,x Qn,x} 



n=2 

A" 



Gaussian -)• - ,_ 1 ^ ct„Tr { i?^ x i? njX } 



n=2 



nonlinear —> + , rc j-i 

^ajTr^}! 2 (Ill) 



^ n=2 



In order to eliminate the nonlinear term in (j 1 1 1 [) , we use the auxiliary variables Q n 

A 



-^£Tr{Q^?} 

n=2 



N K 



Gaussian -—^ J> n Tr {QCgtgCg} 

iV * n=2 

linear +— ^ ]T « n ReTr -Tr {i? n , x } ggtfit I 

* n=2 ^ J 

1 K 

— ^c^Tr^x}! 2 (112) 



cancels out — >■ 



1 n=2 



All nonlinear terms cancel out and are replaced by Gaussian and linear terms. The applicability 
of the efficient updating algorithms discussed in Section [T] follows immediately. 
With auxiliary variables, the partition function of the CYM theory becomes: 



^CYM 



[a, p F ) 



j m [U] im[R,Q) exp (-S F (fa [U]) + ^^ReTr{/tfi x } j (113) 



where /xh[^] is the usual product of S'L r (iV)-invariant Haar measures of the link variables, 
and n[R, Q] is the product of Gaussian measures of the auxiliary variables, 

k / n-l \ 

»[r, q] = n n frniRn*) ^, n (Cx) n ( m ) 

x£A ± n=2 \ m=l / 

with cr's given by: 

al = 2a 2 n = 2a 2 ^ n = a 2 n ^ n = a 2 K ^ = ^- } 2 < n < K, \<m<n-2 (115) 

and / x is a matrix factor that encodes all the information about the deformation terms, which 
is obtained from the linear terms generated by the auxiliary variables: 

h = axsS+E^fgfc^ + E^ (116) 

n=2 \ m=l / 
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The exponent of the Boltzmann factor in the partition function is now fully linear. Hence 
each link variable can be updated with respect to a probability distribution ([6]), where 
is analogous to the typical sum of 'staples'. If the link variable to be updated is parallel to 
a non-compact direction, then V^x contains only the contribution V^ x from the neighboring 
plaquettes; if the link variable is parallel to the compact direction, then it will also contain 
the contribution V^ def from the deformation terms: 

f V*,, if p = l,...,d-l 

r,,, = ( (ii7) 



The plaquette contribution V^ x is given by (26), which already takes into account the possi- 
bility of a fully reduced compact direction (i.e. N t = 1). On the other hand, the contribution 
from the double-trace terms is given by: 

K def = (vU u ^J\ ■f--( v fl u ^ + is] ( 118 ) 

\ 8=0 / \ i=t+l J 

This term resembles the Hermitian conjugate of the Polyakov loop wrapping the compact 
direction and starting at x — (x, t), except that the link variable Ud, x is replaced with f±. 
When implementing the algorithm, one may find convenient to redefine the Polyakov loop f2 x 
to start at x before the link Ud, x is updated, 

N t -l 

^x — > V Y\ ^d,x+i'd = U d x+td ~ ■ ■ ■ Ud,-x+{Nt-l)dP d <*yd,-x.+d ' ' ' ^d,x+(t-l)d (H^) 
i=0 

where i' = i + t (mod Nt). The 'staple' contribution coming from the deformation terms is 
then given by: 

/ m-i \ t 

K dei = /x-bn^j ( 12 °) 

rdei 



If the compact direction is fully reduced, then V x is simply given by: 

K def = / x (121) 



The MC algorithm for lattice CYM theories on R x S is summarised in Appendix A. 3 



4.2 Numerical tests 

We performed numerical simulations of the MC algorithm described above and compared it 
with Metropolis. We considered the lattice regularised version of 577(5) CYM theory defined 
on R 3 x S 1 . We simulated the theory on two different lattices, namely on 10 3 1 and 10 3 3 lattices. 
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The purpose was to test the new algorithm in both situations of an unreduced and a fully 
reduced compact direction. 

For the gauge group SU(5), the CYM action has [_5 / 2J =2 distinct deformation terms. 
The a n parameters attached to them supposedly interpolate between different phases of the 
Yang-Mills theory at fixed coupling. In particular, for large values of fip, the CYM theory 
should interpolate between the deconfining regime for small a n (where the centre symmetry 
associated with the compact direction is spontaneously broken) and a confining regime for 
large a n (where the centre symmetry is intact). This behaviour is a result of the competition 
between the Sp and the S^ef terms in the CYM action. For this reason we chose to perform 
our simulations at (3p = 25.0, which is located in the deconfining (Z^r-broken) regime of the 
fundamental Wilson action. For the 10 3 1 lattice, we chose values for a n that would put the 
CYM theory on two different phases. Specifically, we chose a 'small' a = (0.20, 0.05) and a 
'large' a = (0.60,0.10). For the 10 3 3 lattice, we only chose one value, a = (2.80,0.20). 

For the thermal Metropolis updates, we used a variant of the (1-hit) Cabibbo-Marinari- 
Metropolis (CMM) algorithm described in [19] and adapted to CYM theories. In our CMM 
algorithm, new link proposals are generated using an appropriately tuned Sf action, and are 
then subjected to an accept/reject step with respect to the full CYM action. The acceptance 
rates are tuned to stay within a range of 40-60%. 

For both the Metropolis and the new algorithm, each configuration update consists of 
only one thermal update (no overrelaxation updates were performed). For each configura- 
tion we evaluated the CYM action Scym, which was then used to estimate its expectation 
value (S'cym)- We also evaluated plaquette and Polyakov loop traces. In each simulation we 
performed 398,000 measurements, after discarding the initial 2,000 configurations for equi- 
libration. The simulation parameters and measured observables, together with their naive 
confidence intervals are shown in Table [3] 

Both algorithms agree in terms of the measured values of (S'cym), as can be seen in the first 
and third rows of Table |3j In the second row there is a clear discrepancy in the last significant 
digits, but this is very likely due to the fact that autocorrelations were not taken into account 
in the evaluation of the confidence intervals. From Fig, 14 it is possible to see that r int is very 
large in this example, which means that the confidence intervals on the second row of Table 
[3] are highly underestimated. 

In terms of autocorrelations, the pseudo-heatbath algorithm for the N t = 1 CYM theory 
performs much better than the (optimally tuned) Metropolis algorithm (see Figs, 13 -16). How- 
ever, in the case N t = 3, the pseudo-heatbath algorithm does not show an improvement over 
Metropolis (see Figs 17-18). A possible reason for this behaviour could be an excessive number 
of auxiliary variables in the N t = 3 CYM theory, whose update could easily undermine the 
efficiency of the pseudo-heatbath algorithm. The smaller the number of auxiliary variables, the 
more efficient and faster the algorithm is. In sum, the new updating algorithm for CYM theo- 
ries on R 3 x S 1 are most efficient when the compact direction is fully reduced. Fortunately, this 
is also the most interesting case for a CYM theory, as long as large N volume-independence 
holds nonperturbatively for any volume. 

Finally, Fig, 19 provides a qualitative check for the compatibility of the results from the 
Metropolis and the pseudo-heatbath algorithms. The graphs show the MC histories of the 
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complex trace of the Polyakov loop, for Metropolis and pseudo-heatbath simulations of the 
SU(5) CYM theory on a 10 3 1 lattice. The simulation with a 'small' deformation parameter 
a = (0.20,0.05) results in a deconfining vacuum ((Tr{fi x }) 7^ 0), as expected for a CYM 
theory whose action is dominated by plaquette terms with a large /3p. The simulation with 
'large' deformation parameters a = (0.60, 0.10) results in a confining vacuum ((Tr{fi x }) = 0), 
as expected for a CYM theory whose action is dominated by the centre-stabilising double-trace 
terms. Both situations were correctly captured by the Metropolis and the pseudo-heatbath 
algorithm, which reinforces the validity of the latter. 

5 Discussion 

In this paper we constructed new algorithms for the update of the link variables for two classes 
of pure lattice gauge theories with nonlinear actions. The theories under consideration were 
(i) pure SU(N) lattice gauge theories with a generic mixed Wilson action, and (ii) the lattice 
regularisation of centre-stabilised SU(N) Yang-Mills (CYM) theories defined on R d_1 x S 1 . 

We used a generalisation of the Fabricius-Haan method of auxiliary variables in order to 
construct such algorithms. By adding enough extra degrees of freedom to the lattice gauge 
theory in question, we were able to get rid of the nonlinear terms in its action, and replace 
them with linear terms on the link variables. In this way it is possible to perform pseudo- 
heatbath, overrelaxation or cooling updates on the links of the nonlinear theories, just like in 
the standard SU(N) lattice gauge theory with the fundamental Wilson action. 

As a test for the accuracy of the new algorithms, we evaluated numerically the expectation 
values of some gauge-invariant observables and compared them with a Metropolis evaluation. 
Both quantitative and qualitative results showed a match between the outputs of the new 
pseudo-heatbath and Metropolis algorithms (modulo the autocorrelation correction of the 
confidence levels). 

We also showed numerically that the new algorithms are more efficient than Metropolis. 
Despite the new updates, in general, being slower than Metropolis hits in CPU time (because of 
the large number of auxiliary variables that some lattice theories require) , they also decorrelate 
the gauge configurations very fast. Therefore, when taking autocorrelations into account, the 
new algorithms perform much better than Metropolis. 

The method of auxiliary variables is rather general and may be applied to other lattice the- 
ories with polynomial dependence on the link variables, for example the 3D Georgi-Glashow 
modelj^] the TEK-reduced 2D principal chiral model, etc.. In relation to the CYM theories 
discussed in this paper, the method can also be extended to the more complicated case of 
centre-stabilised Yang-Mills theories defined on manifolds with multiple compactified direc- 
tions, namely ¥L d ~ k x (S r ) k . In particular, it could be used to construct an efficient MC 
algorithm for the zero- volume limit of CYM theories compactified on a d-torus, also known as 
deformed Eguchi-Kawai (DEK) models. An efficient algorithm for the DEK model would be 
very helpful in establishing nonperturbatively (via numerical simulations) if its centre symme- 
try is indeed intact, as claimed by Unsal and Yaffe [13], thus providing the first problem-free 

8 We thank G. Bali for this suggestion. 
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matrix model representation of the planar sector of pure SU(N) Yang-Mills theories on M. d . 
We leave the construction of appropriate algorithms and the study of these theories to later 
publications. 
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A Monte Carlo algorithms 

Here we summarise the MC algorithms for the different SU(N) lattice gauge theories discussed 
in the paper. 

A.l Reduced lattices 

MC algorithm for the fundamental Wilson action ([!]) on partially reduced lattices: 
1. For each reduced plaquette (/if, x): 



(a) Generate one random complex N x N matrix, Qfj,u,x, with normal distribution [i\{Q 



flU,X ) 1 



using the Box-Miiller transform (22). 
(b) Construct new auxiliary variables: 



/ N \ 5 ~ 

Qiw.x — ( 7l ) Quu,x ~\~ Uu %Ul/ x-\-fi ~f~ U v xUu x+0 (1 22) 

\PfJ 

2. For each link variable U^ x '■ 

(a) Construct the sum of 'staples': 

a d 

= ~N E {y^'^^+^l^+p, + Ulx-pU 'n,x-pU v ,x-v+li 

v = \ 

where Q^ >x = Q Uf i,x- 

(b) Update U^ x with respect to the probability distribution ^ using the Cabibbo-Marinari 
pseudo-heatbath algorithm. 

A. 2 Mixed actions 



MC algorithm for generic mixed Wilson actions (28) on lattices with unreduced directions. 



The case of partially reduced directions is discussed in Section |3.5 
1. For each plaquette p = (fiis, x): 

(a) Construct U p = U^ x U UiX+ yU ] ^ x+ Jjt yX 

(b) For each irreducible representation 1Z contributing to the mixed Wilson action: 
i. Generate n-ji random complex N x N matrices Qy with the normal distribution 



jUi(Qp^), using the Box-Miiller transform (22). 



34 



ii. Construct new auxiliary variables using the appropriate /3^, h$ (see Table [lj): 

Q? = -4^Q? + (124) 

iii. Construct f3 using the appropriate gp (see Table 1): 



1=1 



(c) Construct f p : 



f p = 2^/4/ 3 



7v 



(125) 



(126) 



n 



2. For each link variable U^ x '■ 

(a) Construct the sum of 'staples': 

d 



Vfj.,x — ^ ^ f f ^v,xUv,xU ^ > x-\-vU\, 3 ._^_n + U^ x _pfiy^ jX —oU^ ) x—uU l , j x—£>+[i] (127) 



i/=l 



where = = / p . 

(b) Update J7^ )X with respect to the probability distribution ([6| using any of the following 
algorithms: 

• Cabibbo-Marinari pseudo-heatbath 

• SU(2) or SU(N) overrelaxation 

• SU(2) or SU(N) cooling 



A. 3 Double-trace deformations 

MC algorithm for lattice CYM theories on R^ -1 x S 1 , with the possibility of a fully reduced 
'compact' direction. 



1. For each i£ A (<■ 



id-l w cl 



x S 1 ) : 



(a) If Nt = 1, do step (1) from the algorithm A.l 
2. For each xeAi 

(a) Construct Q x = U d ^U d ^ +r ■ ■ U d ^ +{Nt _ 1)S 

(b) Generate a random complex N x N matrix, -Ri )X , with normal distribution H\(R\ x ) 5 using 



the Box-Miiller transform (22) 



35 



(c) Construct the new auxiliary variables: 



Ri* = (siV^J ^i,x + O x --T¥{fi x }l (128) 



(d) Construct = R hx - ±Tr {R hx } 1. 

(e) For n = 2,...,iT : 

i. Construct f2™ 

ii. Generate random complex N x N matrices i?n iX and Q^x with normal distributions 



£*i(-Rn,x) and /ii(Qn^x)) respectively, using the Box-Muller transform (22). 

iii. Construct the new auxiliary variables: 

Rn,* = I^H Rn,* + ^x - ^M^x}l (129) 
1 

= 2 Qg + ^TV^x}^ (130) 

iv. Construct Q^° x = i2 n)X - ^Tr {i? n , x } 1. 

v. For m = 1, . . . , n — 1 : 

A. Generate a random complex N x N matrix with the normal distribution 



Mi(Qrvc), using the Box-Miiller transform (22). 



B. Construct the new auxiliary variables: 

qs = (^^y &+Q^ i) ^+^- m cm) 

(f) Construct / x : 

K ( n-1 1 \ 

/x = aioS + E «M ^ 1} + E oS^Sr 1 ' + 4 Tr {^x} ggt (132) 



n=2 \ m=l 



3. For each link variable Up, x : 



(a) Construct the sum of 'staples' V^ jX as in the step (2a) of the algorithm A.l and multiply 
it by )'/' -V: 



V^ x <- j^rV^ (133) 



(b) If [i = d, then : 
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i. Construct V$ ei : 



/t-1 \T / Nt-l \T 

if N t > i, v; def = p n *W ■/«>!! ^ ( 134 ) 



\ i=0 / \ i=t+l / 

If N t = 1, F x def = / x (135) 

where x = (x, t). 
ii. Add V£ e{ to U x . : 

2iV 

nJ 

(c) Update with respect to the probability distribution ([6| using any of the following 
algorithms: 

• Cabibbo-Marinari pseudo-heatbath 

• SU(2) or SU(N) overrelaxation 

• SU(2) or SU{N) cooling 



K <- V x + ^V x dei (136) 
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Table 1: List of the expressions for /4 and g p 1 ' , useful in the construction of MC algorithms 
for mixed Wilson actions with plaquettes in irreducible representations of SU(N) with N- 
ality k < 3. The symbol '+' stands for 'symmetric representation', either (2) or (3), while 
' — ' stands for 'antisymmetric representation', either (1, 1) or (1, 1, 1). In iV-ality k > 2, two 
different algorithms are suggested for each representation; they differ in the number of auxiliary 
variables. For the adjoint representation, the cases of positive and negative (5a are considered 
separately. Here, p labels plaquettes in the hypercubic lattice, a denotes the sign of the lattice 
coupling, a = sgn(/37e) = ±1, and /3' n is the redefined lattice coupling. 
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Table 2: Expectation values of plaquette characters in the l F + 7V mixed Wilson action, using 
a Cabibbo-Marinari-Metropolis algorithm (M) and a pseudo-heatbath algorithm proposed in 
this paper (H). For the mixed fundamental/adjoint Wilson action, we also include plaquette 
values that were calculated by Hasenbusch and Necco with a Cabibbo-Marinari-Metropolis 
algorithm |19| . The statistical errors in the table do not take autocorrelations into account, 
which is likely the reason for some discrepancies in the last significant digits. 
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Table 3: Expectation values of the SU(5) CYM action for different compactification radii and 
a n , using a Cabibbo-Marinari-Metropolis algorithm (M) and the pseudo-heatbath algorithm 
(H) proposed in this paper. The statistical errors in the table do not take autocorrelations into 
account, which is likely the reason for some discrepancies in the last significant digits. The 
discrepancy in the second row is due to a very slow equilibration of the Metropolis algorithm 
due to large correlations (see Fig. 15). 
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Figure 1: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpv, for the 
SU(3) l F + A' mixed Wilson action, simulated on a 12 3 8 lattice with ((3 F , (3 A ) = (4.00, 2.00), 
via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius-Haan-type pseudo- 
heatbath algorithm (circles). 
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Figure 2: Estimator of the t-dependent integrated autocorrelation time, T int (t), vs. the 
MC time, £mc; f° r the SU(3) l F + A mixed Wilson action, simulated on a 12 3 8 lattice 
with (Pf,Pa) — (4.00,2.00), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 3: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpxj, for 
the SU(4) l F + (2)' mixed Wilson action, simulated on a 8 4 lattice with (Pf,^(2)) — 
(10.665,0.3556), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius- 
Haan-type pseudo- heatbath algorithm (circles). 
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Figure 4: Estimator of the ^-dependent integrated autocorrelation time, r int (t), vs. the MC 
time, £mc; f° r the SU(A) l F + (2)' mixed Wilson action, simulated on a 8 4 lattice with 
{Pf,P(2)) — (10.665,0.3556), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 5: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpu, f° r 
the SU(A) l F + (1,1)' mixed Wilson action, simulated on a 8 4 lattice with (Pf,P(i,i)) = 
(10.665,1.0668), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius- 
Haan-type pseudo- heatbath algorithm (circles). 
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Figure 6: Estimator of the i-dependent integrated autocorrelation time, r int (t), vs. the MC 
time, £mcj ^ or ^ ne SU(A) 'F + (1,1)' mixed Wilson action, simulated on a 8 4 lattice with 
(J3p, /3(h)) = (10.665,1.0668), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 7: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpxj, for 
the SU(5) l F + (3)' mixed Wilson action, simulated on a 8 4 lattice with (/3f,P(3)) = 
(16.665,4.1487), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius- 
Haan-type pseudo-heatbath algorithm (circles). 
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Figure 8: Estimator of the i-dependent integrated autocorrelation time, r int (t), vs. the MC 
time, £mc; f° r the SU(5) l F + (3)' mixed Wilson action, simulated on a 8 4 lattice with 
(/$f,/3(3)) = (16.665,4.1487), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 9: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpv, for the 
SU(5) 'F + (1,1,1)' mixed Wilson action, simulated on a 8 4 lattice with (/3p, /5n 1 1)) = 
(16.665,1.1853), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius- 
Haan-type pseudo- heatbath algorithm (circles). 
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Figure 10: Estimator of the t-dependent integrated autocorrelation time, r; n t(t), vs. the MC 
time, £mco f° r the SU(5) l F + (1, 1, 1)' mixed Wilson action, simulated on a 8 4 lattice with 
(/3p, i = (16.665, 1.1853), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 11: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpu ; for 
the SU(5) l F + (2,1)' mixed Wilson action, simulated on a 8 4 lattice with (Pf,P(2,i)) — 
(16.665,4.7413), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius- 
Haan-type pseudo- heatbath algorithm (circles). 
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Figure 12: Estimator of the t-dependent integrated autocorrelation time, r; n t(t), vs. the MC 
time, £mcj f° r the SU(5) 'F + (2,1)' mixed Wilson action, simulated on a 8 4 lattice with 
(J3p, /3(2,i)) = (16.665,4.7413), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 13: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpu ; for the 
517(5) CYM theory, simulated on a 10 3 1 lattice with (3 F = 25.0 and (a ll a 2 ) = (0.20,0.05), 
via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius-Haan-type pseudo- 
heatbath algorithm (circles). 
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Figure 14: Estimator of the t-dependent integrated autocorrelation time, r int (t), vs. the 
MC time, t MC , for the 5*17(5) CYM theory, simulated on a 10 3 1 lattice with (3 F = 25.0 
and (0:1,02) = (0.20,0.05), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 15: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpu ; for the 
517(5) CYM theory, simulated on a 10 3 1 lattice with (3 F = 25.0 and (a ll a 2 ) = (0.60,0.10), 
via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius-Haan-type pseudo- 
heatbath algorithm (circles). 
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Figure 16: Estimator of the t-dependent integrated autocorrelation time, r int (t), vs. the 
MC time, t MC , for the 517(5) CYM theory, simulated on a 10 3 1 lattice with (3 F = 25.0 
and (ai,Q2) = (0.60,0.10), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 
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Figure 17: Estimator of the autocorrelation function, C(t), vs. the CPU time, tcpu ; for the 
SU(5) CYM theory, simulated on a 10 3 3 lattice with f3 F = 25.0 and (a 1 ,a 2 ) = (2.80,0.20), 
via a Cabibbo-Marinari-Metropolis algorithm (squares) and a Fabricius-Haan-type pseudo- 
heatbath algorithm (circles). 
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Figure 18: Estimator of the t-dependent integrated autocorrelation time, Ti nt (£), vs. the 
MC time, ^mc; f° r the SU(5) CYM theory, simulated on a 10 3 3 lattice with /3p = 25.0 
and (0:1,02) = (2.80,0.20), via a Cabibbo-Marinari-Metropolis algorithm (squares) and a 
Fabricius-Haan-type pseudo-heatbath algorithm (circles). 



50 



a 



0.4 



0.2 



0.0 



-0.2 



-0.4 



S = (0.60,0.10), confining 
S = (0.20, 0.05), deconfining 



confining 



deconfining 



0.0 0.2 0.4 0.6 

iReTr{n x } 



0.8 



P 



0.4 



0.2 



0.0 



-0.2 



-0.4 





a = (0.60, 0.10), confining 






a = (0.20, 0.05), deconfining - 






confining 






5^ ? 












deconfining 











0.0 0.2 0.4 0.6 
iReTr{n x } 



0.8 



Figure 19: Trace of the holonomy, f2 x , around the compact direction of R 3 x S l , in simulations 
of the SU(5) CYM theory on a 10 3 1 lattice with (3p = 25.0, using either a Fabricius-Haan-type 
pseudo-heatbath algorithm (left) or a Cabibbo-Marinari-Metropolis algorithm (right). 
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